<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
<html>
 <head>
  <title>chetrd.f</title>
 <meta name="generator" content="emacs 21.3.1; htmlfontify 0.20">
<style type="text/css"><!-- 
body { background: rgb(255, 255, 255);  color: rgb(0, 0, 0);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: none; }
span.default   { background: rgb(255, 255, 255);  color: rgb(0, 0, 0);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: none; }
span.default a { background: rgb(255, 255, 255);  color: rgb(0, 0, 0);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: underline; }
span.string   { color: rgb(188, 143, 143);  background: rgb(255, 255, 255);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: none; }
span.string a { color: rgb(188, 143, 143);  background: rgb(255, 255, 255);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: underline; }
span.comment   { color: rgb(178, 34, 34);  background: rgb(255, 255, 255);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: none; }
span.comment a { color: rgb(178, 34, 34);  background: rgb(255, 255, 255);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: underline; }
 --></style>

 </head>
  <body>

<pre>
      SUBROUTINE <a name="CHETRD.1"></a><a href="chetrd.f.html#CHETRD.1">CHETRD</a>( UPLO, N, A, LDA, D, E, TAU, WORK, LWORK, INFO )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  -- LAPACK routine (version 3.1) --
</span><span class="comment">*</span><span class="comment">     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
</span><span class="comment">*</span><span class="comment">     November 2006
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     .. Scalar Arguments ..
</span>      CHARACTER          UPLO
      INTEGER            INFO, LDA, LWORK, N
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Array Arguments ..
</span>      REAL               D( * ), E( * )
      COMPLEX            A( LDA, * ), TAU( * ), WORK( * )
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  Purpose
</span><span class="comment">*</span><span class="comment">  =======
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  <a name="CHETRD.19"></a><a href="chetrd.f.html#CHETRD.1">CHETRD</a> reduces a complex Hermitian matrix A to real symmetric
</span><span class="comment">*</span><span class="comment">  tridiagonal form T by a unitary similarity transformation:
</span><span class="comment">*</span><span class="comment">  Q**H * A * Q = T.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  Arguments
</span><span class="comment">*</span><span class="comment">  =========
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  UPLO    (input) CHARACTER*1
</span><span class="comment">*</span><span class="comment">          = 'U':  Upper triangle of A is stored;
</span><span class="comment">*</span><span class="comment">          = 'L':  Lower triangle of A is stored.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  N       (input) INTEGER
</span><span class="comment">*</span><span class="comment">          The order of the matrix A.  N &gt;= 0.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  A       (input/output) COMPLEX array, dimension (LDA,N)
</span><span class="comment">*</span><span class="comment">          On entry, the Hermitian matrix A.  If UPLO = 'U', the leading
</span><span class="comment">*</span><span class="comment">          N-by-N upper triangular part of A contains the upper
</span><span class="comment">*</span><span class="comment">          triangular part of the matrix A, and the strictly lower
</span><span class="comment">*</span><span class="comment">          triangular part of A is not referenced.  If UPLO = 'L', the
</span><span class="comment">*</span><span class="comment">          leading N-by-N lower triangular part of A contains the lower
</span><span class="comment">*</span><span class="comment">          triangular part of the matrix A, and the strictly upper
</span><span class="comment">*</span><span class="comment">          triangular part of A is not referenced.
</span><span class="comment">*</span><span class="comment">          On exit, if UPLO = 'U', the diagonal and first superdiagonal
</span><span class="comment">*</span><span class="comment">          of A are overwritten by the corresponding elements of the
</span><span class="comment">*</span><span class="comment">          tridiagonal matrix T, and the elements above the first
</span><span class="comment">*</span><span class="comment">          superdiagonal, with the array TAU, represent the unitary
</span><span class="comment">*</span><span class="comment">          matrix Q as a product of elementary reflectors; if UPLO
</span><span class="comment">*</span><span class="comment">          = 'L', the diagonal and first subdiagonal of A are over-
</span><span class="comment">*</span><span class="comment">          written by the corresponding elements of the tridiagonal
</span><span class="comment">*</span><span class="comment">          matrix T, and the elements below the first subdiagonal, with
</span><span class="comment">*</span><span class="comment">          the array TAU, represent the unitary matrix Q as a product
</span><span class="comment">*</span><span class="comment">          of elementary reflectors. See Further Details.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  LDA     (input) INTEGER
</span><span class="comment">*</span><span class="comment">          The leading dimension of the array A.  LDA &gt;= max(1,N).
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  D       (output) REAL array, dimension (N)
</span><span class="comment">*</span><span class="comment">          The diagonal elements of the tridiagonal matrix T:
</span><span class="comment">*</span><span class="comment">          D(i) = A(i,i).
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  E       (output) REAL array, dimension (N-1)
</span><span class="comment">*</span><span class="comment">          The off-diagonal elements of the tridiagonal matrix T:
</span><span class="comment">*</span><span class="comment">          E(i) = A(i,i+1) if UPLO = 'U', E(i) = A(i+1,i) if UPLO = 'L'.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  TAU     (output) COMPLEX array, dimension (N-1)
</span><span class="comment">*</span><span class="comment">          The scalar factors of the elementary reflectors (see Further
</span><span class="comment">*</span><span class="comment">          Details).
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  WORK    (workspace/output) COMPLEX array, dimension (MAX(1,LWORK))
</span><span class="comment">*</span><span class="comment">          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  LWORK   (input) INTEGER
</span><span class="comment">*</span><span class="comment">          The dimension of the array WORK.  LWORK &gt;= 1.
</span><span class="comment">*</span><span class="comment">          For optimum performance LWORK &gt;= N*NB, where NB is the
</span><span class="comment">*</span><span class="comment">          optimal blocksize.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">          If LWORK = -1, then a workspace query is assumed; the routine
</span><span class="comment">*</span><span class="comment">          only calculates the optimal size of the WORK array, returns
</span><span class="comment">*</span><span class="comment">          this value as the first entry of the WORK array, and no error
</span><span class="comment">*</span><span class="comment">          message related to LWORK is issued by <a name="XERBLA.78"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  INFO    (output) INTEGER
</span><span class="comment">*</span><span class="comment">          = 0:  successful exit
</span><span class="comment">*</span><span class="comment">          &lt; 0:  if INFO = -i, the i-th argument had an illegal value
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  Further Details
</span><span class="comment">*</span><span class="comment">  ===============
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  If UPLO = 'U', the matrix Q is represented as a product of elementary
</span><span class="comment">*</span><span class="comment">  reflectors
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Q = H(n-1) . . . H(2) H(1).
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  Each H(i) has the form
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     H(i) = I - tau * v * v'
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  where tau is a complex scalar, and v is a complex vector with
</span><span class="comment">*</span><span class="comment">  v(i+1:n) = 0 and v(i) = 1; v(1:i-1) is stored on exit in
</span><span class="comment">*</span><span class="comment">  A(1:i-1,i+1), and tau in TAU(i).
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  If UPLO = 'L', the matrix Q is represented as a product of elementary
</span><span class="comment">*</span><span class="comment">  reflectors
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Q = H(1) H(2) . . . H(n-1).
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  Each H(i) has the form
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     H(i) = I - tau * v * v'
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  where tau is a complex scalar, and v is a complex vector with
</span><span class="comment">*</span><span class="comment">  v(1:i) = 0 and v(i+1) = 1; v(i+2:n) is stored on exit in A(i+2:n,i),
</span><span class="comment">*</span><span class="comment">  and tau in TAU(i).
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  The contents of A on exit are illustrated by the following examples
</span><span class="comment">*</span><span class="comment">  with n = 5:
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  if UPLO = 'U':                       if UPLO = 'L':
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">    (  d   e   v2  v3  v4 )              (  d                  )
</span><span class="comment">*</span><span class="comment">    (      d   e   v3  v4 )              (  e   d              )
</span><span class="comment">*</span><span class="comment">    (          d   e   v4 )              (  v1  e   d          )
</span><span class="comment">*</span><span class="comment">    (              d   e  )              (  v1  v2  e   d      )
</span><span class="comment">*</span><span class="comment">    (                  d  )              (  v1  v2  v3  e   d  )
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  where d and e denote diagonal and off-diagonal elements of T, and vi
</span><span class="comment">*</span><span class="comment">  denotes an element of the vector defining H(i).
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  =====================================================================
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     .. Parameters ..
</span>      REAL               ONE
      PARAMETER          ( ONE = 1.0E+0 )
      COMPLEX            CONE
      PARAMETER          ( CONE = ( 1.0E+0, 0.0E+0 ) )
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Local Scalars ..
</span>      LOGICAL            LQUERY, UPPER
      INTEGER            I, IINFO, IWS, J, KK, LDWORK, LWKOPT, NB,
     $                   NBMIN, NX
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Subroutines ..
</span>      EXTERNAL           CHER2K, <a name="CHETD2.141"></a><a href="chetd2.f.html#CHETD2.1">CHETD2</a>, <a name="CLATRD.141"></a><a href="clatrd.f.html#CLATRD.1">CLATRD</a>, <a name="XERBLA.141"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Intrinsic Functions ..
</span>      INTRINSIC          MAX
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Functions ..
</span>      LOGICAL            <a name="LSAME.147"></a><a href="lsame.f.html#LSAME.1">LSAME</a>
      INTEGER            <a name="ILAENV.148"></a><a href="ilaenv.f.html#ILAENV.1">ILAENV</a>
      EXTERNAL           <a name="LSAME.149"></a><a href="lsame.f.html#LSAME.1">LSAME</a>, <a name="ILAENV.149"></a><a href="ilaenv.f.html#ILAENV.1">ILAENV</a>
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Executable Statements ..
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Test the input parameters
</span><span class="comment">*</span><span class="comment">
</span>      INFO = 0
      UPPER = <a name="LSAME.156"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( UPLO, <span class="string">'U'</span> )
      LQUERY = ( LWORK.EQ.-1 )
      IF( .NOT.UPPER .AND. .NOT.<a name="LSAME.158"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( UPLO, <span class="string">'L'</span> ) ) THEN
         INFO = -1
      ELSE IF( N.LT.0 ) THEN
         INFO = -2
      ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
         INFO = -4
      ELSE IF( LWORK.LT.1 .AND. .NOT.LQUERY ) THEN
         INFO = -9
      END IF
<span class="comment">*</span><span class="comment">
</span>      IF( INFO.EQ.0 ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Determine the block size.
</span><span class="comment">*</span><span class="comment">
</span>         NB = <a name="ILAENV.172"></a><a href="ilaenv.f.html#ILAENV.1">ILAENV</a>( 1, <span class="string">'<a name="CHETRD.172"></a><a href="chetrd.f.html#CHETRD.1">CHETRD</a>'</span>, UPLO, N, -1, -1, -1 )
         LWKOPT = N*NB
         WORK( 1 ) = LWKOPT
      END IF
<span class="comment">*</span><span class="comment">
</span>      IF( INFO.NE.0 ) THEN
         CALL <a name="XERBLA.178"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>( <span class="string">'<a name="CHETRD.178"></a><a href="chetrd.f.html#CHETRD.1">CHETRD</a>'</span>, -INFO )
         RETURN
      ELSE IF( LQUERY ) THEN
         RETURN
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Quick return if possible
</span><span class="comment">*</span><span class="comment">
</span>      IF( N.EQ.0 ) THEN
         WORK( 1 ) = 1
         RETURN
      END IF
<span class="comment">*</span><span class="comment">
</span>      NX = N
      IWS = 1
      IF( NB.GT.1 .AND. NB.LT.N ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Determine when to cross over from blocked to unblocked code
</span><span class="comment">*</span><span class="comment">        (last block is always handled by unblocked code).
</span><span class="comment">*</span><span class="comment">
</span>         NX = MAX( NB, <a name="ILAENV.198"></a><a href="ilaenv.f.html#ILAENV.1">ILAENV</a>( 3, <span class="string">'<a name="CHETRD.198"></a><a href="chetrd.f.html#CHETRD.1">CHETRD</a>'</span>, UPLO, N, -1, -1, -1 ) )
         IF( NX.LT.N ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Determine if workspace is large enough for blocked code.
</span><span class="comment">*</span><span class="comment">
</span>            LDWORK = N
            IWS = LDWORK*NB
            IF( LWORK.LT.IWS ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Not enough workspace to use optimal NB:  determine the
</span><span class="comment">*</span><span class="comment">              minimum value of NB, and reduce NB or force use of
</span><span class="comment">*</span><span class="comment">              unblocked code by setting NX = N.
</span><span class="comment">*</span><span class="comment">
</span>               NB = MAX( LWORK / LDWORK, 1 )
               NBMIN = <a name="ILAENV.212"></a><a href="ilaenv.f.html#ILAENV.1">ILAENV</a>( 2, <span class="string">'<a name="CHETRD.212"></a><a href="chetrd.f.html#CHETRD.1">CHETRD</a>'</span>, UPLO, N, -1, -1, -1 )
               IF( NB.LT.NBMIN )
     $            NX = N
            END IF
         ELSE
            NX = N
         END IF
      ELSE
         NB = 1
      END IF
<span class="comment">*</span><span class="comment">
</span>      IF( UPPER ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Reduce the upper triangle of A.
</span><span class="comment">*</span><span class="comment">        Columns 1:kk are handled by the unblocked method.
</span><span class="comment">*</span><span class="comment">
</span>         KK = N - ( ( N-NX+NB-1 ) / NB )*NB
         DO 20 I = N - NB + 1, KK + 1, -NB
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Reduce columns i:i+nb-1 to tridiagonal form and form the
</span><span class="comment">*</span><span class="comment">           matrix W which is needed to update the unreduced part of
</span><span class="comment">*</span><span class="comment">           the matrix
</span><span class="comment">*</span><span class="comment">
</span>            CALL <a name="CLATRD.235"></a><a href="clatrd.f.html#CLATRD.1">CLATRD</a>( UPLO, I+NB-1, NB, A, LDA, E, TAU, WORK,
     $                   LDWORK )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Update the unreduced submatrix A(1:i-1,1:i-1), using an
</span><span class="comment">*</span><span class="comment">           update of the form:  A := A - V*W' - W*V'
</span><span class="comment">*</span><span class="comment">
</span>            CALL CHER2K( UPLO, <span class="string">'No transpose'</span>, I-1, NB, -CONE,
     $                   A( 1, I ), LDA, WORK, LDWORK, ONE, A, LDA )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Copy superdiagonal elements back into A, and diagonal
</span><span class="comment">*</span><span class="comment">           elements into D
</span><span class="comment">*</span><span class="comment">
</span>            DO 10 J = I, I + NB - 1
               A( J-1, J ) = E( J-1 )
               D( J ) = A( J, J )
   10       CONTINUE
   20    CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Use unblocked code to reduce the last or only block
</span><span class="comment">*</span><span class="comment">
</span>         CALL <a name="CHETD2.255"></a><a href="chetd2.f.html#CHETD2.1">CHETD2</a>( UPLO, KK, A, LDA, D, E, TAU, IINFO )
      ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Reduce the lower triangle of A
</span><span class="comment">*</span><span class="comment">
</span>         DO 40 I = 1, N - NX, NB
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Reduce columns i:i+nb-1 to tridiagonal form and form the
</span><span class="comment">*</span><span class="comment">           matrix W which is needed to update the unreduced part of
</span><span class="comment">*</span><span class="comment">           the matrix
</span><span class="comment">*</span><span class="comment">
</span>            CALL <a name="CLATRD.266"></a><a href="clatrd.f.html#CLATRD.1">CLATRD</a>( UPLO, N-I+1, NB, A( I, I ), LDA, E( I ),
     $                   TAU( I ), WORK, LDWORK )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Update the unreduced submatrix A(i+nb:n,i+nb:n), using
</span><span class="comment">*</span><span class="comment">           an update of the form:  A := A - V*W' - W*V'
</span><span class="comment">*</span><span class="comment">
</span>            CALL CHER2K( UPLO, <span class="string">'No transpose'</span>, N-I-NB+1, NB, -CONE,
     $                   A( I+NB, I ), LDA, WORK( NB+1 ), LDWORK, ONE,
     $                   A( I+NB, I+NB ), LDA )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Copy subdiagonal elements back into A, and diagonal
</span><span class="comment">*</span><span class="comment">           elements into D
</span><span class="comment">*</span><span class="comment">
</span>            DO 30 J = I, I + NB - 1
               A( J+1, J ) = E( J )
               D( J ) = A( J, J )
   30       CONTINUE
   40    CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Use unblocked code to reduce the last or only block
</span><span class="comment">*</span><span class="comment">
</span>         CALL <a name="CHETD2.287"></a><a href="chetd2.f.html#CHETD2.1">CHETD2</a>( UPLO, N-I+1, A( I, I ), LDA, D( I ), E( I ),
     $                TAU( I ), IINFO )
      END IF
<span class="comment">*</span><span class="comment">
</span>      WORK( 1 ) = LWKOPT
      RETURN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     End of <a name="CHETRD.294"></a><a href="chetrd.f.html#CHETRD.1">CHETRD</a>
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>
